
## PEASANT RESISTANCE AND ECONOMIC AFFLUENCE: LESSONS FROM PARAGUAY
## Main regression tables and plots (paper)
## January 2024

## Code Authors:
## Germán Feierherd (gfeierherd@udesa.edu.ar) 
## Jorge Mangonnet (jorge.mangonnet@vanderbilt.edu)


#### Preamble ------------------------------------------------------------------

## Workspace setup
rm(list = ls())		         # clear all objects in memory
options(max.print = 5000)  # expand display
options(scipen = 999)      # disable sci notation
dev.off()                  # reload graphic device
cat("\014")                # clear console

## Load/install packages
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  readxl,
  interflex,
  marginaleffects,
  haven,
  readr,
  reshape2,
  magrittr,
  dplyr,
  plm,
  lfe,
  pcse,
  lmtest,
  stargazer,
  xtable,
  texreg,
  estimatr, 
  fixest,
  margins, 
  sjPlot, 
  multiwayvcov,
#  rgdal,
  ggplot2,
  wesanderson, 
  fixest,
  did,
  sf)

## Set working directory
gf <- "~/Dropbox/Working_Papers/investigacion protesta Paraguay/"; main <- gf  # german's path
jm <- "~/Dropbox/Research/investigacion protesta Paraguay/"; main <- jm        # jorge's path
setwd(main)

## Functions
source(paste0(main, "4_code/functions_paraguay.R"))

## Set CRS for shapefiles 
wgs84 <- CRS("+proj=longlat +datum=WGS84 +no_defs")

## Load full dataset
load(paste0(main, "4_code/DataSetPY.Rda"))


#### Regression table ----------------------------------------------------------

#### Table 1: Peasant Resistance to Land Encroachment, 2000-2014 ----
#### Note: 'coeftest' for cluster standard errors

# Model 1: Prices x Suitability, w/o controls
m1 <- plm(ln_conflict ~ I(ln_suitability*ln_price) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m1) 
se1 <- coeftest(m1, vcov = vcovHC(m1, type = "HC2", cluster = "group")); se1

# Model 2: Prices x Suitability, w/ controls
m2 <- plm(ln_conflict ~ I(ln_suitability*ln_price) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m2) 
se2 <- coeftest(m2, vcov = vcovHC(m2, type = "HC2", cluster = "group")); se2

# Model 3: Prices x Suitability x Settlements, w/o controls
m3 <- plm(ln_conflict ~ I(ln_suitability*ln_price) + I(ln_price*subsistence_farms_median) + I(ln_suitability*ln_price*subsistence_farms_median) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m3) 
se3 <- coeftest(m3, vcov = vcovHC(m3, type = "HC2", cluster = "group")); se3

# Model 4: Prices x Suitability x Settlements, w/ controls
m4 <- plm(ln_conflict ~ I(ln_suitability*ln_price) + I(ln_price*subsistence_farms_median) + I(ln_suitability*ln_price*subsistence_farms_median) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m4) 
se4 <- coeftest(m4, vcov = vcovHC(m4, type = "HC2", cluster = "group")); se4

# Model 5: Prices x Suitability x Committees, w/o controls
m5 <- plm(ln_conflict ~ I(ln_suitability*ln_price) + I(ln_price*orgs_locales_median) + I(ln_suitability*ln_price*orgs_locales_median) + factor(year), index = c("codigo"), model = "within", effect = "individual", na.action = na.omit, data = dta); summary(m5) 
se5 <- coeftest(m5, vcov = vcovHC(m5, type = "HC2", cluster = "group")); se5

# Model 6: Prices x Suitability x Committees, w/ controls
m6 <- plm(ln_conflict ~ I(ln_suitability*ln_price) + I(ln_price*orgs_locales_median) + I(ln_suitability*ln_price*orgs_locales_median) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", na.action = na.omit, data = dta); summary(m6) 
se6 <- coeftest(m6, vcov = vcovHC(m6, type = "HC2", cluster = "group")); se6

# Create table
stargazer(
  m1, m2, m3, m4, m5, m6, 
  se = list(se1[,2], se2[,2], se3[,2], se4[,2], se5[,2], se6[,2]), 
  p = list(se1[,4], se2[,4], se3[,4], se4[,4], se5[,4], se6[,4]),
  omit = c("factor", "Constant", "distance_hubs", "untitled_share", "ln_pop"),
  omit.stat = c("adj.rsq", "ser", "f"),
  covariate.labels = c("Prices$_{t-1}$ $\\times$ Suitability", 
                       "Prices$_{t-1}$ $\\times$ Settlements", 
                       "Prices$_{t-1}$ $\\times$ Suitability $\\times$ Settlements", 
                       "Prices$_{t-1}$ $\\times$ Committees", 
                       "Prices$_{t-1}$ $\\times$ Suitability $\\times$ Committees"),
  add.lines = list(c("","","","","",""), 
                   c("Controls", "No", "Yes", "No", "Yes", "No", "Yes")), 
  dep.var.labels.include = FALSE,
  model.names = FALSE,
  df = FALSE, 
  multicolumn = TRUE, 
  type = "latex",
  star.cutoffs = c(0.1, 0.05, 0.01),
  star.char = c("*", "**", "***"),
  float = TRUE, 
  float.env = "sidewaystable",
  style = "ajps", 
  title = "Peasant Resistance to Land Encroachment, 2000-2014", 
  label = "main",
  notes.append = FALSE,
  notes = c("\\parbox[t]{20cm}{\\textit{Note}: All models include municipality and year fixed effects. Standard errors clustered by municipality in parentheses.}",
            "$^{*}p<.1$, $^{**}p<.05$, $^{***}p<.01$"),
  out = paste0(main,"1_writing/2023/tables/table_main.tex")
)

# Remove
remove(m1, m2, m3, m4, m5, m6, se1, se2, se3, se4, se5, se6)


#### Marginal effects plots ----------------------------------------------------

#### Figure 3: Marginal Effect of Agricultural Prices on Peasant Resistance by Land Suitability ----

# Based on model 1 (Table 1)
m1 <- plm(ln_conflict ~ ln_suitability*ln_price + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta)

# Covariance matrix for cluster standard errors
covMat1 <- vcovHC(m1, type = "HC2", cluster = "group")

# Plot
pdf(file = paste0(main, "1_writing/2023/figures/inter_main.pdf"), width = 5, height = 5)
interaction_plot_continuous(m1, effect="ln_price", moderator="ln_suitability", interaction="ln_suitability:ln_price", varcov=covMat1, minimum=0.75, maximum=1.4, incr="default", num_points = 50, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=T, title=NULL, xlabel="Land Suitability", ylabel="Marginal Effect of Agricultural Prices")
dev.off()


#### Figure 4: Marginal Effect of Agricultural Prices on Peasant Resistance by Land Suitability, Subsistence Settlements, and Peasant Committees ----

# Based on model 3 & 5 (Table 1)
m3 <- plm(ln_conflict ~ ln_suitability*ln_price + ln_price*subsistence_farms_median + ln_suitability*ln_price*subsistence_farms_median + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta)
m5 <- plm(ln_conflict ~ ln_suitability*ln_price + ln_price*orgs_locales_median + ln_suitability*ln_price*orgs_locales_median + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta)

# Covariance matrices for cluster standard errors
covMat3 <- vcovHC(m3, type = "HC2", cluster = "group")
covMat5 <- vcovHC(m5, type = "HC2", cluster = "group")

# Plot (a): subsistence settlements
pdf(file = paste0(main, "1_writing/2023/figures/inter_subsist.pdf"), width = 5, height = 5)
interaction_plot_threeway_ylim(m3, effect="ln_price", moderator1="ln_suitability", moderator2="subsistence_farms_median", varcov=covMat3, minimum=0.75, maximum=1.4, incr="default", num_points = 100, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=F, title=NULL, xlabel="Land Suitability", ylabel="Marginal Effect of Agricultural Prices")
dev.off()
# Plot (b): peasant committees
pdf(file = paste0(main, "1_writing/2023/figures/inter_organiz.pdf"), width = 5, height = 5)
interaction_plot_threeway_ylim(m5, effect="ln_price", moderator1="ln_suitability", moderator2="orgs_locales_median", varcov=covMat5, minimum=0.75, maximum=1.4, incr="default", num_points = 50, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=F, title=NULL, xlabel="Land Suitability", ylabel="Marginal Effect of Agriculutral Prices")
dev.off()


#### Figure 5: Conditional Adjusted Predictions ----
#### Note: based on model 3 & 5 (Table 1)

# Visual options
theme_set(theme_minimal())
colores <- c(adjustcolor("blue", alpha.f = 0.5), adjustcolor("red", alpha.f = 0.5))
options(ggplot2.discrete.fill = colores)
options(ggplot2.discrete.colour = colores)
options(width = 1000)


## (a) Subsistence settlements --

# Fit model (just for illustrating predicted values)
fit_subsist <- lm(ln_conflict ~ ln_suitability * ln_price * subsistence_farms_median + factor(year), data = dta)

# Plot predictions
cp_subsist <- 
  plot_predictions(fit_subsist,
  # ses
  vcov = vcovHC(fit_subsist, type = "HC2", cluster = "group"), 
  # set predictor values 
  condition = list(
    "ln_suitability" = "minmax",
    "subsistence_farms_median" = unique,
    "ln_price" = "threenum"),
  draw = TRUE) +
  ## ggplot configuration ##
  # theme
  theme_bw() + 
  theme(axis.text = element_text(size = 15, hjust = 1),
        axis.title = element_text(size = 15, hjust = .5),
        strip.text = element_text(size = 13), 
        legend.text = element_text(size = 13)) +
  # labels
  labs(x = "Suitability", y = "Events of Peasant Resistance (IHS)") + 
  coord_cartesian(ylim = c(-0.1, 1.5)) + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "black", size = .45) + 
  # legends & titles
  scale_colour_discrete(name = "", labels = c("Low","High")) + 
  guides(color = guide_legend(override.aes = list(linetype = 0, size = 1.5))) + 
  facet_grid(. ~ ln_price, labeller = as_labeller(c(`-SD` = "Low Prices: -SD", `Mean` = "Mean Prices", `+SD` = "High Prices: +SD")))

pdf(file = paste0(main, "1_writing/2023/figures/condp_subsist.pdf"))
cp_subsist
dev.off()


## (b) Peasant committees --

# Fit model (just for illustrating predicted values)
fit_orgs <- lm(ln_conflict ~ ln_suitability * ln_price * orgs_locales_median + factor(year), data = dta)

# Plot predictions
cp_orgs <- 
  plot_predictions(fit_orgs,
  # ses
  vcov = vcovHC(fit_orgs, type = "HC2", cluster = "group"), 
  # set predictor values 
  condition = list(
    "ln_suitability" = "minmax",
    "orgs_locales_median" = unique,
    "ln_price" = "threenum"),
     draw = TRUE) +
  ## ggplot configuration ##
  # theme
  theme_bw() + 
  theme(axis.text = element_text(size = 15, hjust = 1),
        axis.title = element_text(size = 15, hjust = .5), 
        strip.text = element_text(size = 13), 
        legend.text = element_text(size = 13)) +
  # labels
  labs(x = "Suitability", y = "Events of Peasant Resistance (IHS)") + 
  coord_cartesian(ylim = c(-0.1, 1.5)) + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "black", size = .45) + 
  # legends & titles
  scale_colour_discrete(name = "", labels = c("Low","High")) + 
  guides(color = guide_legend(override.aes = list(linetype = 0, size = 1.5))) + 
  facet_grid(. ~ ln_price, labeller = as_labeller(c(`-SD` = "Low Prices: -SD", `Mean` = "Mean Prices", `+SD` = "High Prices: +SD")))

pdf(file = paste0(main, "1_writing/2023/figures/condp_organiz.pdf"))
cp_orgs
dev.off()


#### Descriptive plots ---------------------------------------------------------

#### Figure 1: Agricultural Prices, Land Suitability, and Peasant Resistance (2000-2014) ----

# Divide munis into higher suit and lower suit
dta$land_suit_bin <- ifelse(dta$land_suitability > median(dta$land_suitability), "Higher Suitability", "Lower Suitability")
# Data frame for plots
d1 <- dta %>% 
  dplyr::select(land_suit_bin,year,n_conflict) %>% 
  dplyr::group_by(land_suit_bin,year) %>% 
  dplyr::summarise(n_conflict = sum(n_conflict))

pdf(file = paste0(main, "1_writing/2023/figures/ts_resist_new.pdf"), width = 5, height = 5)
barplot(d1$n_conflict ~ d1$land_suit_bin + d1$year, ylim = c(0,250), col = c(adjustcolor("red", alpha.f = 0.50), adjustcolor("blue", alpha.f = 0.50)),
        xlab = "Year", ylab = "Events of Peasant Resistance \nAgricultural Prices  (USD per metric ton)", border = "white", 
        legend.text = unique(d1$land_suit_bin),
        args.legend = list(x = "topright", bty = "n"))
dev.off()

pdf(file = paste0(main, "1_writing/2023/figures/ts_prices_new.pdf"),  width = 5, height = 5)
with(dta %>% dplyr::select(year,commodity_price) %>% distinct(), plot(commodity_price ~ year, type="l", ylab="", axes=F))
axis(side = 1)
axis(side = 4, seq(0, 250, 50))
mtext("", side = 4, line = 3)
dev.off()

pdf(file = paste0(main, "1_writing/2023/figures/ts_prices_resist_new.pdf"),  width = 6, height = 5)

par(mar = c(4, 6, 4, 4) + .1)  # Leave space for z axis

barplot(d1$n_conflict ~ d1$land_suit_bin + d1$year, ylim = c(0,250), col = c(adjustcolor("red", alpha.f = 0.50), adjustcolor("blue", alpha.f = 0.50)),
        xlab = "Year", ylab = "Events of Peasant Resistance", border="white", 
        legend.text = unique(d1$land_suit_bin),
        args.legend = list(x = "topleft", bty = "n"))

axis(side = 4, seq(0, 250, 50))
mtext("Agricultural Prices  (USD per metric ton)", side = 4, line = 3)

par(new = TRUE)

with(dta %>% dplyr::select(year, commodity_price) %>% distinct(), plot(x = year, y = commodity_price, type = "l", ylab = "", xlab = "", axes = FALSE))

dev.off()


#### Figure 2: Peasant Resistance and Land Suitability by Municipality ----

# Group the data
cs <- dta %>%
  dplyr::select(codigo, n_conflict, ln_suitability, population) %>% 
  group_by(codigo) %>%
  summarise(n_conflict = sum(n_conflict),
            ln_suitability = mean(ln_suitability),
            population = min(population)) %>%
  mutate(ln_conflict = inv_hyp(n_conflict), 
         pc_conflict = (n_conflict / population) * 10000)

# Load municipal polygons and merge
py <- st_read(paste0(main, "7_gis/paraguay_2002_municipios"), "paraguay_2002_municipios") %>%
  st_transform(wgs84) %>%
  merge(y = cs, by = "codigo") %>%
  mutate(id = seq(1, 223, 1))

# Country border
py_border <- st_boundary(st_union(py))

## (a) Peasant resistance --
map.resist <- 
  ggplot() +
  theme_void() + 
  geom_sf(data = py, aes(fill = ln_conflict), color = "transparent", size = 1) +
  geom_sf(data = py_border, color = "black", fill = NA, size = 1) +
  scale_fill_gradient(name = "Events (IHS)", low = 'white', high = "red", na.value = "white") + 
  theme(legend.position = "right") + 
  theme(legend.text = element_text(size = 16)) +
  theme(legend.title = element_text(size = 16)) + 
  theme(plot.title = element_text(size = 18)) + 
  theme(plot.title = element_text(hjust = 0.5))

png(file = paste0(main, "1_writing/2023/figures/map_resist.png"), width = 8, height = 8, res = 300, units='in')
map.resist
dev.off()

## (b) Land suitability --
map.suit <- 
  ggplot() +
  theme_void() + 
  geom_sf(data = py, aes(fill = ln_suitability), color = "transparent", size = 1) +
  geom_sf(data = py_border, color = "black", fill = NA, size = 1) +
  scale_fill_gradient(name = "Tons per\nhectare (log)", low = 'white', high = "blue", na.value = "white") + 
  theme(legend.position = "right") + 
  theme(legend.text = element_text(size = 16)) +
  theme(legend.title = element_text(size = 16)) + 
  theme(plot.title = element_text(size = 18)) + 
  theme(plot.title = element_text(hjust = 0.5))

png(file = paste0(main, "1_writing/2023/figures/map_suit.png"), width = 8, height = 8, res = 300, units='in')
map.suit
dev.off()
